#include "prs_sum.h"

double_t prs_sum(const float_t* sumData, const size_t sumLength)
{
	double_t sum = 0.0f;
	for (int64_t i = 0; i < (int64_t)sumLength; i++)
	{
		sum += (double_t)logf(sqrtf(sumData[i]));
	}
	return sum;
}

double_t prs_sum_omp(const float_t* sumData, const size_t sumLength)
{
	double_t sum = 0.0f;
	#pragma omp parallel for reduction(+:sum)
	for (int64_t i = 0; i < (int64_t)sumLength; i++)
	{
		sum += (double_t)logf(sqrtf(sumData[i]));
	}
	return sum;
}

double_t prs_sum_avx(const float_t* sumData, const size_t sumLength)
{
	#define BlockWidth 8    // 块宽. AVX寄存器能一次处理8个float.

	double_t sum = 0.0f;

	size_t blockCnt = sumLength / BlockWidth;		// 块数.
	size_t remainCnt = sumLength % BlockWidth;		// 块数.
	__m256 avxSum = _mm256_setzero_ps();			// 求和变量。[AVX] 赋初值0
	

	const float_t* p = sumData,*q= (const float_t*)&avxSum;						// AVX批量处理时所用的指针.
	for (int64_t i = 0; i < (int64_t)blockCnt; i++)
	{
		avxSum = _mm256_add_ps(avxSum, _mm256_log_ps(_mm256_sqrt_ps(_mm256_load_ps(p))));    // [AVX] 单精浮点紧缩加法
		
		p += BlockWidth;
	}
	for (size_t i = 0; i < BlockWidth; i++)
	{
		sum += (double_t)q[i];
	}
	for (size_t i = 0; i < remainCnt; i++)
	{
		sum += p[i];
	}

	return (double_t)sum;
}

double_t prs_sum_omp_avx(const float_t* sumData, const size_t sumLength)
{
	#define procCnt 1600
	double_t subSum[procCnt];
	double_t sum=0.0;
	#pragma omp parallel for
	for (int32_t i = 0; i < procCnt; i++)
	{
		subSum[i] = prs_sum_avx(sumData + sumLength * i / procCnt, sumLength / procCnt);
	}
	for (size_t i = 0; i < procCnt; i++)
	{
		sum += subSum[i];
	}
	return sum;
}

double_t prs_sum_merge(const double_t op1, const double_t op2)
{
	return op1 + op2;
}